Seismic data can be used to interpret or to infer sub-surface geology, making it useful for the location, identification and exploitation of petroleum and minerals. However since seismic traces are often contaminated by random noise, seismic data sets typically undergo a series of conventional statistical processes (known as “seismic processing”) before the data can be so used. It is advantageous to remove noise at an early stage in processing, since this improves the ability to perform all subsequent processing work. Conventional seismic processing disadvantageously degrades or distorts signal prior to removing sufficient noise for some purposes. Conventional seismic processing methods are based on functions and their transforms that it is often useful to think of as occupying two domains. These domains have been referred to as the function and transform domains, but more commonly they are known as the time and frequency domains. Operations performed in one domain have corresponding operations in the other. For example, the convolution operation in the time domain becomes a multiplication operation in the frequency domain and the reverse is also true, permitting users to move easily between domains so that operations can be performed where they are easiest. A seismic audio signal being sampled at 8 Hz means that at each successive eighth of a second a measurement of the intensity of the signal is taken. The Fourier transform decomposes or separates a waveform or function into sinusoids of different frequency, which sinusoids sum to the original waveform. It identifies or distinguishes the different frequency sinusoids and their respective amplitudes. The Discrete Fourier Transform (DFT) is useful because a digital computer works only with discrete data, numerical computation of the Fourier transform of f(t) requires discrete sample values of f(t), which are called fk. In addition, a computer can compute the transform F(s) only at discrete values of s, in other words it can only provide discrete samples of the transform, Fr. It is well understood that a Discrete Fourier Transform (“DFT”) is a Fourier transform calculated for a wavelet over a finite interval so that values are given only for the fundamental frequency (reciprocal of the interval) and its harmonics. And a “trace” is a time domain record of seismic data from a single electromagnetic channel.
Within the seismic industry sections of data are comprised of several “grids” that form cubes of data when used to correlate Cartesian (x-y surface) coordinates to either a time domain or frequency domain that covers a finite range of time or frequency respectively. In a Time Cube; a trace is a temporal record of data from a single seismic channel, the sensor for which is located on a point of intersection defined according to a given scale of abscissas (i.e. along the X-axis) and ordinates (i.e. along the Y-axis)—each of which points of intersection are identified by a co-ordinate pair and all of which points of intersection form a planar rectangular “grid”. The 3rd dimension being time, a “trace” forms as seismic energy (typically in the form of acoustic signal amplitude) sample recordings are captured from the sensor location identified by the subject co-ordinate pair, over a finite range of time (e.g. 2 seconds) at defined (e.g. 10 millisecond) intervals. The scope and depth of the data comprising a Time Cube is related to the number of elements (Cartesian co-ordinate pairs), the time, and the sampling rate through that time. For example, in each 25×25 “grid” there will be 625 traces forming a time record of signals (comprising both genuine reflected energy and noise) recorded in relation to the subject section (i.e. 1 data element associated with each co-ordinate pair=625 data elements per sample period). With each trace having samples captured every 10 milliseconds over 2 seconds, this Time Cube would contain 125,000 data elements.
Prior Art FIG. 1 illustrates a typical geological exploration setup for acquiring seismic data. Positioned at or just below the earth's surface 110, an energy source 120 (typically explosive high-energy or vibrational low energy source having an effective length of only 20 to 40 ms) generates at least one sound wave or signal (the wavefront resulting from which is typically recorded for approx. 3 seconds) having sufficient energy to follow path 130 down into the earth a suitable distance to reflect at the interface of any changes in geology (commonly known as “events” or “reflectors”) 140, the reflected energy traveling back to the surface via path 150 is simultaneously recorded by receivers 160 (commonly geophones positioned as an array for 3D, or a line for 2D exploration). In marine seismic the sound wave is generated just below the surface of the water and the reflected energy detected by hydrophones. For each such sound wave generation, or “shot”, the reflected signal returning to the surface via path 150 creates a “trace”, which is a single recording at a receiver 160. Traces are detected and recorded in the form of a time series of sample measurements of particle velocity (for land data) or pressure change (for marine data). Many shots are taken to generate a seismic data set, often resulting in hundreds of millions of traces that may be stacked of summed in a variety of ways. When high energy impulsive seismic sources are used, the level of the detected true earth response seismic signal is usually greater than the ambient noise. However, when low energy surface seismic sources are used, the ambient noise can be at a level greater than the true earth response seismic signal. For this reason, seismic-trace recordings are often made involving the repeated initiation of a low energy surface seismic source at about the same origination point, thereby producing a sequence of seismic-trace data based on seismic wave reflections and/or refractions that have traveled over about the same path and therefore have approximately the same travel times. The process of adding such seismic-trace data together for improving the signal-to-noise ratio of the composite seismic-trace recording is known as “vertical compositing” or “vertical stacking.” It should be distinguished from “horizontal stacking,” a process applied to a sequence of seismic-trace data based on seismic wave reflections from approximately the same subsurface point (referred to as the “common-depth point,” or “CMP”) but which has been generated and recorded at different surface locations. Horizontal stacking of CMP seismic-trace data requires that time corrections (called “normal moveout,” or “NMO,” corrections) be applied before the traces are summed together, since travel times from seismic source to detector are unequal for each trace in the sequence. It can be assumed that the true earth response seismic signal embedded in each trace is coherent and in phase (correlated) and that the noise is random and incoherent (uncorrelated) with zero mean value. In general, the objective of vertical stacking is to maximize the signal-to-noise ratio of the resultant recording. Reflectors that are not “flat” are said to “dip” or slope.
The use of a low energy vibrator can be more economical than the use of dynamite. Furthermore, as compared to the use of a high-energy impulsive seismic source, such as dynamite, the frequency of the seismic waves generated by a vibrator can be selected by controlling the frequency of the pilot signal to the power source, such as a hydraulic motor, which drives the vibrator. More particularly, the frequency of the pilot signal to the vibrator power source can be varied, that is, “swept,” for obtaining seismic-trace data at different frequencies. A low energy seismic wave, such as generated by a vibrator, can be used effectively for seismic prospecting if the frequency of the vibrator “chirp” signal which generates the seismic wave is swept according to a known pilot signal and the detected seismic wave reflections and/or refractions are cross correlated with the pilot signal in order to produce seismic-trace recordings similar to those which would have been produced with a high energy impulsive seismic source. Typically, the pilot signal is a swept frequency sine wave that causes the vibrator power source to drive the vibrator for coupling a swept sine wave “chirp” signal into the earth. The swept frequency operation yields seismic-trace data that enables different earth responses to be analyzed, providing a basis on which to define the structure, such as the depth and thickness, of the subsurface formations. It is a problem that recorded seismic-trace data always includes some background (ambient) noise in addition to the detected seismic waves reflected and/or refracted from the subsurface formations (referred to as the “true earth response”). Noise can be classified as “stationary” and “non-stationary”, both of which can be random. Stationary noise is random noise such as atmospheric electromagnetic disturbances that are statistically time-invariant over the period of acquisition of seismic-trace data for producing a recording. Non-stationary noise is random and often occurs as bursts or spikes generally caused by wind, traffic, recorder electrical noise, et cetera, which are statistically time-variant over the period of acquisition of seismic-trace data for producing a recording and exhibits relatively large excursions in amplitude. In connection with swept frequency operation of low energy vibrator seismic prospecting, it is common practice to vertically stack, or sum, the seismic-trace data from a series of initiations, that is, sequential swept frequency operations, to produce a composite seismic-trace recording for the purpose of improving the signal-to-noise ratio of the seismic-trace data. Unfortunately, the commonly used technique of vertically stacking trace data is inadequate in the presence of non-stationary noise that appears during such seismic prospecting.
Seismic data is acquired in two principal geometries: 2-D and 3-D. In 2-D acquisition, shots and receivers are positioned along a (not necessarily straight) surface line. In 3-D acquisition, shots and receivers are positioned over a 2-D surface area. Seismic data related to 3D geologic volumes necessarily includes random noise that may be isolated from the signal data to different degrees by different conventional techniques, including an eigenimage filtering technique that is 2D in nature and disadvantageously does not account for additional available information respecting the formation.
For 2-D acquisition the main product of seismic processing is a 2-D stacked “section” (illustrated in Prior Art FIG. 2), one of the dimensions representing horizontal position along acquisition line 210, and the other dimension representing time 220. For 3-D acquisition, seismic processing resulting in a 3-D stacked section (illustrated in Prior Art FIG. 3), two of the representing edges 310 and 320 of the acquisition surface area, and the other dimension representing time 330.
Known seismic processing steps include common-midpoint (CMP) stacking, where traces are collected into groups having roughly the same midpoints between the locations of the shot by which they were generated and the receiver at which they were detected. For each recorded time sample, the magnitudes or values of every trace in the group are summed together, producing a single “stacked” trace for each group. Such stacking commonly reduces the amount of data that must be processed by a factor of 10 to 100.
Geological interpretation is easiest and most successful on seismic sections having low levels of noise, and thus one of the objects of seismic processing is to remove as much noise as possible. Noise can be broadly categorised as random, coherent, or monochromatic. Random noise may be defined as noise that is uncorrelated between traces and spectrally broad band. Some of the causes of random noise are the effects of wind and other disruptions on the seismic receiver and cable, poor penetration of seismic energy through the earth (particularly in the near surface beneath the shot or receiver), and numerous natural and manmade seismic energy sources apart from the intended one. The most common stage to carry out random-noise removal is after CMP stacking. A number of methods have been developed to do this, including: f-K transform (March and Bailey, 1983), f-x prediction (Canales, 1984; Soubaras, 1994), Karhunen-Loeve transform (Jones and Levy, 1987; Al-Yahya, 1991), eigenimage (Ulrych, Sacchi, and Freire, 1999), spectral matrix filtering (Gounon, Marse, and Goncalves, 1998), and Radon transform (Russell, Hampson, Chun, 1990(a) and 1990(b)).
The foregoing methods work on 2-D data sets, but can be adapted for 3-D stacked sections by slicing the data volume along one of be spatial dimensions, filtering each of these slices separately as if it were a 2-D section, and then recomposing the 3-D volume. This can then be repeated in the opposite spatial direction. Such methods are not optimum in that they fail to fully exit the large amount of data available within a short radius of each spatial point of the 3-D volume. For this reason, “true 3-D” methods have been developed that work in both spatial dimensions simultaneously, including: f-xy prediction (Chase, 1992; Soubaras, 2000) and f-kk transform (Peardon and Bacon, 1992).
Random noise removal before CMP stacking is less common. There are, however, at least two advantages to removing random noise as early in the processing stream as possible. First, it improves the performance of subsequent processes, notably deconvolution, statics correction, and velocity analysis. Second, it has the potential to be more effective since more data is available before stacking, providing better statistical redundancy. At the same time, extra data means that random noise removal before CMP stacking requires more computation. Noise removal before statics or deconvolution faces the problem of “surface consistent effects”, meaning effects that are constant within each shot and receiver, but that may change radically even between adjacent shots and receivers. If these effects have not been corrected before noise removal then the noise removal process must preserve them, one method for this is surface-consistent f-x prediction (Wang, 1996).
Another application of noise removal is common-offset or common-angle stacks for amplitude versus offset (AVO) or amplitude versus angle (AVA) analysis. Such stacks are used for the automatic computation of parameters for interpretation. These stacks require a low level of noise so the computed parameters are as accurate as possible. In 2-D acquisition, AVO/AVA stacks form a 3-D volume in which the two spatial dimensions are CMP and either offset or angle. In the offset-angle dimension there may be only one or two dozen traces, such that much of the data is on or near a spatial boundary. Disadvantageously even the better noise removal methods, such as f-xy prediction, do not perform well near spatial boundaries, resulting in distortion of the signal, and possible distortion of the computed parameters—creating the need for a noise removal method that performs well at spatial boundaries.
Some conventional noise removal methods such as Karhunen-Loeve, time-domain eigenimage filtering, Radon transform, and f-K and f-kk transforms, only work well on plains-type data received from geological formations in which most of the reflectors are flat. Disadvantageously on more structured data received from geological formations in which reflectors are strongly dipping or sloping, these methods become computationally expensive, difficult to use, or less effective. Accordingly, there is a need for a method of noise removal that performs well on all types of geology.
Many methods for removing noise from seismic data are implemented using matrix operations. For example, matrix compression is one process for determining an approximation to a given matrix, which approximation consumes less space than the original matrix. However, matrix rank reduction is to determine the nearest (with respect to a particular matrix norm) rank-K matrix to a given matrix. A matrix norm measures the size of a matrix, for example, the “Frobenius norm” is the square root of the sum of the square of the matrix elements, whereas the “L1 matrix norm” is the sum of the absolute values of the elements of the given matrix.
Prior art reviewed includes U.S. Pat. No. 5,379,268 to Hutson, who teaches compression of the subject matrix (see claim 1(b)) as the core of his improvement Compression and rank reduction have some similarities, but are not identical, in that rank reduction can be used as a step within matrix compression, but rank reduction can also be performed without any resulting matrix compression. Hutson in 268 teaches compression—by active decomposition (expanding the original data matrix), then actively and selectively zeroing (whereas modifying to a value other than zero does not “compress”) singular values in one of the resulting matrices, thereafter recomposing those matrices into a single matrix that approximates the original—after which 268 proceeds to teach (see claim 1(c)(1)) further processing of the signals approximated by the compressed matrix. However, 268 is vague even about the kind of compression and processing performed. For example modifying without zeroing cannot compress while zeroing inappropriate selections can actually be counterproductive by eliminating signal rather than noise.
Disadvantageously, existing noise removal algorithms do not handle erratic noise well. For example, both SVD and Lanczos methods of matrix rank reduction attempt to find a rank-K matrix that is nearest to the input matrix in a Frobenius-norm sense. This is appropriate for removing random noise that has a Gaussian, or bell-shaped, statistical distribution—which does not perform as well when erratic non-Gaussian noise bursts are present.
Processing seismic data is typically time consuming and expensive because it involves large quantities of data. Therefore, it is desirable to provide a solution to at least some of the above-described problems of the prior art reducing either or both the quantity of data processed or the amount of processing required in relation to that data. The prior art includes some matrix rank reduction and frequency domain based methods, but there is a need for a new combination of such methods that overcomes the above disadvantages, particularly for the purpose of noise suppression in sections of seismic data.